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Abstract. We present a comprehensive analytical model of aeolian sand transport 
in saltation. It quantifies the momentum transfer from the wind to the transported 
sand by providing expressions for the thickness of the saltation layer and the apparent 
surface roughness. These expressions are for the first time entirely derived from basic 
physical principles. The model further predicts the sand transport rate (mass flux) 
and the impact threshold shear velocity. We show that the model predictions are in 
very good agreement with experiments and numerical state of the art simulations of 
aeolian saltation. 
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1. Introduction 

Saltation is the dominant mechanism of aeohan sand transport on Earth's deserts under 
turbulent wind flow. Unidirectional wind accelerates sand grains, which perform hops 
of typical shapes. During their hops the wind continuously transfers momentum to the 
grains. Therefore the wind momentum decreases, resulting in reduced wind velocities 
not only within, but also above the saltation layer. Above the saltation layer the average 
horizontal wind velocity profile u{z) follows the well known law [1], 



where is the wind shear velocity, k = 0.4 is the von Karman constant, and z* is the 
apparent roughness of the moving saltation layer. In the absence of sand transport z* 
becomes Zo, which is the surface roughness of a quiescent sand bed. In the presence of 
sand transport the magnitude of z* depends on how much momentum is absorbed by 
the saltation layer. It is crucial for the development of sand transport models, but also 
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for landscape modelers and coastal managers to know z* as function of and Zg- For 
instance in aeolian dune models z* is a key quantity in the computation of the wind 
field over a non-fiat topography, in which the shear velocity varies with the spatial 
position [2, 3, 4, 5]. The main purpose of this paper is to derive a novel prediction of 
zl, which is entirely based on physical principles. 

Deriving a scaling law for was also approached by previous studies. First Owen 
[6] suggested 

z: oc ^, (2) 

where g is the gravity constant. (2) is also known as the Charnock relation, since 
Charnock [7] derived (2) for the roughness of a wind-blown water surface. Owen 
[6] based his formula on the assumptions that the average lift-off velocity vi, with 
which a grain leaves the bed, is proportional to m*, that opposing drag forces can be 
neglected, and that z* is proportional to the average hop height of grains h, also called 
saltation height. However the author's assumptions fail to agree with measurements. 
Experimental studies [8, 9, 10, 11] found Vi and consequently h to be almost independent 
of within the measured range. Furthermore Zq cannot be proportional to h, since 
was found to strongly vary with in experiments [12, 13, 14] in contrast to h. In 
addition Sherman [15] found that (2) leads to strong discrepancies with experiments 
close to the impact threshold Ut, which is the threshold shear velocity at which saltation 
can be sustained through the splash process. On Earth Ut is below the fluid entrainment 
shear velocity, needed to entrain sand grains from the soil by fluid lift [1, 16]. Sherman 
[15] therefore extended (2) to the so-called modified Charnock relation, 

i^^i^. (3) 

9 

which ensures that at the threshold = Ut, where no particles are moving, the 
roughness is unchanged, z* = Zq. Although (3) was successfully validated with the data 
set of Sherman and Farrell [14] for z*, it shares the same lack of physical justification 
as (2). 

A much more physical approach was presented by Raupach [17]. From the mixing 
length approximation [18], the author derived 

ln^= fl-^Vn^^, (4) 
Zo \ u^J l.YSZo 

where Zg is the decay height of the grain shear stress profile Tg{z), which the author 
assumed to be exponentially decreasing, 

Tg{z) = Tgoe-'/'% (5) 

where Tgo = 'Tg{0). Tg{z) describes how much momentum is transferred, at each height z, 
from the fluid to the grains per unit soil area and time. The difficulty in the usage of (4) 
is the undetermined quantity Zg. Raupach [17] therefore assumed Zs to be proportional 
to the saltation height h, which he in turn assumed to be proportional to ul/g like before 
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Owen [6]. However, as we discussed before, Owen's assumption is in disagreement with 
experiments. Relations very similar to (4) were also obtained in two further studies 
[19, 20]. These studies achieved agreement with the experimental data of Rasmussen et 
al. [12], by introducing an ad-hoc fit relation for Zg- Andreotti [19] found that the data 
set can be well fitted if Zg scales with \/d, where d is the mean particle diameter, and 
he therefore suggested 

Zs (X ^/sgd ( , (6) 

where s — Ps/ Pw is the ratio of sand density Ps and fluid density p^, and // is the 
kinematic viscosity. A scaling law very similar to (6) was also used by Duran and 
Herrmann [20]. However (6) is very weakly founded on physics. Its only justification 
is the resulting agreement of (4) with the data set of Rasmussen et al. [12]. Therefore 
there is a great necessity to either validate (6) or to derive a new expression for Zg from 
physical principles. 

Within this study wc do the latter. Wc present a comprehensive analytical model 
of aeolian saltation, which aims to significantly improve previous analytical models 
[1, 6, 21, 22, 23]. It for the first time provides expressions for Zg and z* entirely derived 
from physical principles. Our analysis will reveal that Zg is a measure for the thickness 
of the saltation layer and not proportional to the average hop height h. The model 
furthermore incorporates expressions for other important sand transport quantities, such 
as the mass flux Q and the impact threshold Ut. The model is based on the concept 
of mean motion, meaning that average quantities are used for its description. In our 
model we separately consider the horizontal and vertical transport of grains. For each 
we analytically balance the average force and work rate per unit soil area applied by the 
wind on a grain during a trajectory versus the respective amounts applied by the soil on 
the grains during an impact. This results in a model parameter a, describing the ratio 
between the average vertical and horizontal force per unit soil area, and another model 
parameter ^, describing the ratio between the average work rate per unit soil area in the 
vertical motion and the horizontal motion. With theoretical arguments it is shown that 
a and ^ are nearly independent of and the atmospheric conditions and only slightly 
varying with the buoyancy-reduced gravity g and the particle diameter d. The model 
further contains a third parameter 7, defined as the ratio between Zg and the effective 
height of the mean motion Zm- The final relations for Zg, z*, and Q are functions of a, (3, 
7, and Uf. We afterwards extend our model in such a way that also Ut can be computed. 
Thereby a fourth parameter 77 comes into play, describing the ratio between the average 
particle velocity, reduced by the particle slip velocity, and the average wind velocity. 
In our model we use the assumption that the grain shear stress profile is exponentially 
decaying, equivalent to what was assumed in previous studies [17, 20] (see (5)). 

We validate our model with the apparent roughness data of Rasmussen et al. [12] 
for five different grain sizes, with combined impact threshold data of several studies 
[24, 25, 26], and with mass ffux data of Creyssels et al. [10]. Furthermore, through 
simulations with the numerical state of the art model of Kok and Renno [27] we 
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support our derived expressions and our statement that the model parameters are nearly 
independent of u^, the atmospheric conditions, as well as g and d. 

The manuscript is structured in the following way. It starts with a comprehensive 
model description in Section 2, which is followed by the model validation in Section 3 
and a discussion of the results in Section 4. The appendices incorporate long calculations 
and side information. There is also a glossary at the end of the manuscript, which helps 
to keep track of the mathematical symbols. 

2. Model description 

It is the main purpose of our paper to derive a novel expression for the apparent 
roughness 2;* during aeolian sand transport in steady state. The main focus of our 
model lies therefore in the analytical description of the momentum and energy transfer 
from the wind to the grains. Momentum and energy transfer are the main causes 
for the increase of the surface roughness Zg of a quiescent sand bed to the apparent 
roughness of a moving saltation layer. In detail we use Newton's law to obtain 
equations, which balance the average force and work rate per unit soil area applied 
during a grains trajectory with the respective amounts applied during an impact. The 
ratio of the average force (work rate) per unit soil area for the vertical motion and force 
(work rate) per unit soil area for the horizontal motion is the definition of our model 
parameter a {(3). After applying the balance laws we show that the decay height Zg 
of the grain shear stress profile Tg{z) and subsequently z^ as well as the mass flux Q 
can be calculated from the impact threshold Ut and our model parameters. Afterwards 
the model is extended in such a way that also Uf can be calculated as function of the 
model parameters. As previous studies [17, 20] we assume an exponentially decreasing 
Tg{z) (see (5)) and extensively use it in our calculations. This assumption is therefore 
discussed in a separate paragraph. 

For the balance laws, we only consider wind drag and gravity as driving forces, but 
neglect turbulent lift forces, the Magnus force, electrostatic forces, and momentum as 
well as energy changes by mid-air collisions between grains for simplicity reasons and 
because gravity and drag dominate the sand transport [27]. Furthermore we simplify 
the description, by only considering average quantities, which implies that we neglect 
turbulent fluctuations of the wind velocities. Further simpliflcations are the use of 
monodisperse, spherical sand grains, being transported above a horizontal sand bed. The 
probably most crucial of all these simpliflcations is the negligence of mid-air collisions 
between saltating grains. The effect of such collisions has only rarely been subject of 
scientiflc studies [28, 29, 30, 31], because for a comparison between saltation with and 
without mid-air collisions, one has to turn off mid-air collisions, what is possible (and 
common) in numerical simulations, but impossible in experiments. According to the 
most recent numerical study (Figure 5 in Ren and Huang [31]), the change of the mass 
flux due to mid-air collisions is less than 10% for 3.5ut. This is below the typical 
measurement error of mass flux measurements (> 10%). According to this study the 
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neglegance of mid-air collisions and therefore our model simplifications are acceptable 
up to at least u^, ~ 3.5ut. 

This section is separated in several subsections. It starts with the presentation of 
notations and definitions, which are used for the description of our model, in Section 2.1. 
In Section 2.2 follows a short discussion of our main model assumption, the exponentially 
decreasing grain shear stress profile. After that the balance laws are appfied, first for 
the force per unit soil area in Section 2.3 and then for the work rate per unit soil area in 
Section 2.4. Subsequently we discuss the invariance of the model parameters a and P in 
Section 2.5. Afterwards we obtain a novel relation for Zg in Section 2.6, which is further 
discussed in Section 2.7. Then in Section 2.8 relations for and Q as a function of Ut 
and the model parameters are obtained. In Section 2.9 the model is extended, in order 
to allow for the computation of Ut as well. 



2.1. Notations and definitions 

For the coming analytical calculations, we henceforth use the following notations: An 
index x refers to the horizontal component of a given quantity, which coincides with the 
direction of the wind, an index z to the vertical direction, whereby z is also the height 
above the sand bed. Furthermore we differentiate between the upward and downward 
part of a grain's trajectory by indices t and I, respectively. Quantities evaluated at the 
sand bed z — incorporate an additional index a. In particular quantities, which refer 
to a grain's impact, consist of the indices a and I, if the quantity is evaluated before 
the impact, and the indices a and t, if the quantity is evaluated after the impact. 

In order to keep the manuscript simple, it is advantageous to predefine quantities, 
which are used in the following calculations. One quantity is the average particle mass 
per unit volume p{z), transported at height z. p{z) integrated over the whole saltation 
layer describes the mass M of transported sand per unit soil area. 

oo 

M = y p{z)dz. (7) 



Since we differentiate between upward and downward movement p{z) can be divided in 
the mass of upward and downward moving particles per unit volume 

piz)^p^{z)+p^iz). (8) 

Other important quantities are the average vectorial wind velocity profile, u{z), whose 
2;-component is zero u{z) :— Ux{z) — \u{z)\, and the average vectorial particle velocity 
profile for the upward (downward) part of the trajectory v^(^i,){z). The difference between 
both velocities is denoted as 

Based on these definitions, we further define the following velocity differences by 

Avx{z) = Vxi{z) - fxt(^), (10) 
Avz{z) = Vzi{z) - Vz^{z), (11) 
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^vUz)^vl{z)-vl^{z), (12) 

^vl(z)-vli(z)-vlr(^), (13) 
as well as the local vertical mass flux (f){z) by 

(j){z) = p^{z)v,^{z) = -pi{z)v,^{z), (14) 

where we used that the vertical upward flux must exactly compensate the downward 
flux in steady state. Note that Vz],{z) and thus Avz{z) are negative. Using (8), (14) can 
be rewritten as 

With these definitions the average gain of horizontal and vertical momentum of a 
transported grain per unit soil area and time between the two times it crosses height z 
can be written as 

Tg(z) = (l)(z)Av^(z) (16) 
for the horizontal and 

Pg{z) = <P{z)Av,{z) (17) 

for the vertical momentum gain per unit soil area and time. Tg{z) is also known as the 
grain shear stress profile [17, 20, 23] and Pg{z) can be seen as a grain normal stress 
(grain pressure) profile. Note that, by inserting (15) in (17), one obtains 

Pg{z) = p{z)v,i{z)v,^{z), (18) 

which is identical to the definition of the granular pressure in previous studies [10, 32], 
if vertical drag is neglected. 



2.2. Grain shear stress profile 

In this section we discuss and motivate our main model assumption of an exponentially 
decreasing grain shear stress profile (see also (5)) 

Tgiz) = ve-^/^^ (19) 

(19) is justified in the following manner. First, an approximately exponentially 
decreasing mass density profile p{z) was measured in wind tunnels [10, 11]. Although 
not necessarily identical, the profiles Tg{z) and p{z) should at least behave in a similar 
manner. Therefore it is very reasonable that also Tg{z) decreases approximately 
exponentially. Second, Tg{z) has been obtained from numerical simulations [27, 33], 
which indeed showed an approximately exponential decrease. This is shown in Figure 
1 for simulation results with the numerical model of Kok and Renno [27]. It should be 
noted that the mass density profile p{z) strongly deviates from the exponential shape 
at very small heights in the simulations. Such a deviation is also present for the grain 
shear stress profile Tg{z) (see Figure 1 at heights very close to zero), however to a much 
lesser extent. 
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Figure 1. Plot of the grain sliear stress profile Tg{z) obtained from numerical 
simulations with the model of Kok and Renno [27] for = 0.25m/ s (blue), = 
0.5m/s (red), — 0.75m/s (green), and = Im/s (brown). The simulations are 
performed under Earth conditions with a mean diameter d = 250/im. Over a large 
part Tg{z) decays exponentially for all shear velocities. 



2.3. Force balance 

As already pointed out, the description of the momentum transfer from the wind to the 
grains is a key ingredient towards a description of the feedback effect of sand transport 
on the wind profile and thus a first step towards a prediction of the apparent roughness 
z*. Newton's second law for grains moving in a particular trajectory, indicated by a 
lower index T, can be written as {dz = viz^(^i)dt) 

dfixt(i) 



Pita)^i^t(;) 



dz 
d^i.t(4) 
dz 



(20) 
(21) 



for the upward (downward) part of this trajectory, where /ixt(4.) fiztil) 
horizontal and vertical components of the total average force fit(4,) per unit volume 
acting on the grain in the upward (downward) part of this trajectory. For the single- 
trajectory case pit^izt = —PuVizi = 0(0) is constant with height z [34]. Summing the 
upward and downward part of (20) and (21), respectively, followed by averaging over 
all trajectories and integration over height therefore yields 



(z)AvJz) 



fA^')dz', 



(22) 



Pgiz) = (l)iz)Av,iz) = / f,iz')dz' 



(23) 
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where we approximated the trajectory average of products as the product of trajectory 
averages in (pAv^ and (f>Avz- Here Av^o = Av^i^O) and Av^o = Avz{0), and fx{z) and 
f^{z) are the averages of /i^ = fi^^ + /i^; and /i^ = fu^ + f^i over all trajectories, 
respectively. The terms on the left hand side of (22) and (23) describe the average 
horizontal and vertical force per unit soil area applied by the soil on the grains during 
an impact and the right hand side the average horizontal and vertical force per unit soil 
area applied by the wind on the grains during a trajectory, respectively. In the next 
steps we evaluate the integrals in (22) and (23). Therefore we first need an expression 
for the total trajectory-averaged force f . Since we neglect the Magnus force, turbulent 
lift forces, and momentum transfer through collisions, f is only composed of the drag 
and the gravity force. Further approximating the trajectory average of products as the 
product of trajectory averages, f can be written as 
3 

f = ^ (PtC<i(M)MVrt + PiCd{Vrl)VriVri) - pQB^, (24) 

where is the unit vector in 2;-direction, g = ^g with s = Ps/pw is the buoyancy- 
reduced gravity (for most atmospheres g ^ g), and Vrt(4.) = |vrt(4.)|- Cd is the drag 
coefficient, which is a function of the particle Reynolds number and therefore of t'rt(4.)- 
For many drag laws in the literature e.g. [35, 36] the dependency of Cd on a velocity 
difference V can be described by a law of the type 

where Co and Coo = Cd{oo) are dimensionless parameters. Such a drag law strongly 
simplifies f^, which becomes 

where (14), (16), and Vrz^d) = —Vz\{i) were used. If a drag law of another type than 
(25) was used e.g. [37], (26) would still be valid in very good approximation. On the 
other hand we rewrite f^ as 

fx = -^{CdMVrVrx) , (27) 

where () denotes a weighted average of a quantity / between the upward and downward 
movement, (/) = {p^f^ + Pif])l P- Now we can evaluate the integral in (22) using the 
expression we derived for f^- We obtain 



r.(0) = / ^^{Cd{v^)v,v^Mz - '^'^^^f-'^ , (28) 



where the overbar and the capital letters denote the average of a quantity / over height, 
F — pfdz/ pdz. We further used (7), the approximation Vr ~ Vrx, which is 
reasonable since in aeolian saltation the horizontal motion dominates the vertical one, 
and we approximated the height average of the products by the products of the height 
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averages {Cd{vr)vrVrx) ^ Cd{Vr)V^. On the other hand Pgo can now be calculated from 
our expression for f^. It becomes 

Pgo = pM = -9M-^-^j^, (29) 

where we used our assumption (19). With the evaluation of the integrals, we can now 
define a' and our first model parameter a as 

/ Pgo ^Vzo /oA\ 

'go '-^^xo 



a — a ^ J — — — — — — zz^i W-l-j 



, SCpo^g _ ^sgd 
^sd ~ ^Ca{Vr)vl 

where we used (16) and (17) as well as the notations /^v^o — At'a;(0) and Av^o = ^Vz{Q). 
The advantage of this definition of a lies in the fact that a' is almost independent of 
li*, atmospheric conditions, g, and d as we show later in a separate chapter. For many 
conditions, we can approximate 

a ^ a', (32) 

since a' is typically much larger than CooZg/ (sd), as will be verified later. It mainly 
means that the gravity force is large in comparison to the vertical drag force, fz ~ —gp- 
We can thus formulate relevant sand transport quantities as a function of a constant a. 
For instance from (31) we obtain a direct relation between the average velocity difference 
Vr and a, writing 

CdiVrWl = (33) 

and further, using (29), (30), and (32), a direct relation between the grain shear stress 
Tgo at the bed and the mass of transported grains per unit soil area M, writing 

Tgo = a-^gM. (34) 



2.4- Work rate balance 

The second important ingredient towards a description of the feedback effect of the 
grain motion on the wind profile and towards a prediction of z* is the description of 
the energy transfer from the fluid to the grains. Since we discuss a purely Newtonian 
problem, we separate the horizontal and vertical motion. The work rate balance with 
respect to the horizontal (vertical) motion can be obtained by multiplying (20) ((21)) 
with fiitU) ("^121(4-))' summing the upward and downward part, integrating over height, 
and averaging over all trajectories. It yields 

oo 

Im^vl = J ifxtvx^ + Uvxi.)dz, (35) 



oo 

Im^vl = J iUvzt + Uvzi)dz, (36) 
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where Av^^ = Ai>^(0), Av"^^ = Af^(O), and we approximated the trajectory average of 
products as the product of trajectory averages in 0(O)Av^p and (f){0)Av'^^. The terms 
on the left hand side of (35) and (36) describe the average work rate during an impact 
and the right hand side the average work rate during a trajectory for the horizontal and 
vertical motion, respectively. Analogous to (24) and (27) we can now write 

/a;tM + fxi^xi = ^ {CdiVr)VrVrxVx) , (37) 

fz^V,^ + f,iv,i = -j^iCd^Vrv"^,), (38) 
where we used (14) in (38). Analogous to (28), integration now approximately yields 
l<f>mvl ^ ^^CdiVrWlv, (39) 

l<t>mvl ^ - ^CaiVrWrV!, (40) 

where V — V x describes the average particle velocity. Note that (40) would write = 0, 
if vertical drag is neglected {Av\^ — and {C^iyr^VrV^^ — 0). Evaluating the integrals, 
we now define ^' as 




(41) 



(41) means that the average granular temperature is proportional to V rV ■ This is 
different from Creyssels et al. [10] who found to be approximately equal to (i;^)(0). 
The main reason for this difference is that the authors neglected vertical drag, whereas 
our description considers it (see (40)). As before for a', the advantage of (41) hes in 
the fact that /3' is almost independent of it*, atmospheric conditions, ^, and d as we will 
show in the following. 



2. 5. Invariance of a' and (3' 

Since our model relations, including the final relation for z^, will be expressed as 
functions of a' and f3', it is important to discuss, how these parameters change with 
varying conditions. As can be seen from (30) and (41), both parameters are ratios 
of certain velocity differences evaluated at the soil z = and therefore related to the 
splash-entrainment process. The splash-entrainment process dominates the entrainment 
of bed grains in aeolian steady state saltation, because the entrainment by wind is small 
due to a strong reduction of the wind velocity close to the sand bed, which even leads 
to decreasing near-surface velocities with increasing [20]. Since fluid-entrainment is 
not relevant, each impacting grain must exactly lead to one grain leaving the surface 
(rebound or ejection of new grains) on average for steady state sand transport. The 
average number grains leaving the surface per impacting grain can however only depend 
on the average impact velocity Vi, angle 9i, and the relevant bed properties g and d. 
Furthermore, for given values of g and d, the average velocity Vi and angle 6i of a grain 
leaving the surface after an impact, called lift-off velocity and angle, can only depend 
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(a) (b) 




0.0001 0.001 0.01 0.0001 0.001 0.01 

gd / (m^/s^) gd / (m^/s^) 



Figure 2. a' (blue) and (red) computed with tlie model of Kok [40] (see Appendix 
A) plotted versus gd. Here the gravity was fixed to either the Earth (solid lines) or 
the Mars value (dashed lines) and d varied. 



on Vi and 6i. This means, for given values of g, d, and 9i, there are unique values Vi, 
Vi, and 6i, which fulfill that one impacting grain makes one grain leave the surface on 
average. 

Another necessary condition, which must hold, is that Vzoi = Vi sin 6i must be 
smaller than Vzo-\ = vi sin 6i due to friction with the air, it would be equal in the absence 
of drag. The validity of this condition was observed in saltation experiments [38]. 
Further, from collision experiments Oger et al. [39] found that this condition is only 
fulfilled for small impact angles 6i ^ 15°. The authors also found that the number 
of ejected particles significantly decreases with decreasing 6i, meaning that the range 
of ^j- values, in which saltation can be sustained, should be rather narrow. In fact, it 
has been measured in experiments [8] that the average impact angle is approximately 
constant, 9i ~ 11°, with increasing u^:, and thus the other splash quantities as well stay 
approximately constant with varying u^. Here and henceforth we refer to situations 
above the impact threshold > Ut and within our model limits (not too large m*), 
when mentioning dependencies on m^,. This means that in particular the parameters a' 
and (3' are approximately independent of u^, and atmospheric properties, and thus only 
functions of g and d. 

In order to get an idea of how a' and f3' behave as functions of g and d, we use 
the model of Kok [40], with which Vzoi and Vzot and thus a' and (3' can be computed. 
The model is briefly described in Appendix A. The model results in functions a'{gd) 
and (3'{gd) plotted in Figure 2. For the plots the gravities of Earth, g = 9.81m/s^, and 
Mars, g = 3.71m/ s"^, were fixed and d varied. As can be seen, even as functions of g 
and d the variance of a' and /3' is small according to this model. Furthermore a' is 
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of the order of unity, what justifies the approximation (32) for many conditions. For 
instance for Earth conditions the neglected term can be estimated as being of the order 
of 1/25, since Zs is the same order of magnitude as the saltation height h, which has 
been measured as being equal to about iOd [11], Coo ~ 1) and s ~ 1000 on Earth. 



2.6. A novel relation for Zg 

The definitions of the parameters a and obtained from the momentum and energy 
balances can now be used to express the decay height z^ of the grain shear stress profile 
Tg{z) as a function of a and (3'. As already pointed out Zs is the key quantity towards 
a prediction of the apparent roughness z^ . For the calculation of Zg we use 

(v!) = ei^f±£i^ = (42) 

Pt + Pi 

where we inserted (14). We further approximate the arithmetic average of \vz\\ and 
by their geometric average 

-^r~ ^ n V l^^tl = ^/-VztVzl, (43) 



2 

what is reasonable since |i'2;-|-|/|i'24.| < l^zot 1/^204. 1 a-nd It'^ot 1/^2:04. 1 is about 1.5 as 
measurements indicate [39], which means that the error of this approximation is less 
than 3%. Using (14), (16), (42) and (43), we express Tg as 

= lp^/{^Av,, (44) 
This allows us to rewrite (22) as 



where we used (19). Since (45) is valid for all heights z, it must be particularly valid 
for the average over height. We therefore approximately calculate Zg as 



^Cd{Vr)VrVr^) ?>Cd{Vr)Vl 9 



where /3 = (3'AVx/{2V) and we used (31) and (41). /3 is our second model parameter and 
like f3' approximately constant, since we expect that Av^ is in leading order proportional 
to V. (46) is the main contribution of our paper, since it is, to our knowledge, the first 
physically based prediction of and therefore the most important part towards a novel 
physically based prediction of z^. If the values of the model parameters a and /3 are 
known, V remains the only undetermined quantity in (46), since Vr can be calculated 
by (33) as a function of a. There is evidence that Zg does not only describe the decay 
of Tg{z), but also the decay of p{z) for large z. This can be seen from (45), which says 
that Tg{z) decays in the same way as p{Cd{vr)VrVrx) ■ Since {Cd{vr)vrVrx) is only slightly 
decaying with height for large z, the decaying behaviors of p{z) and Tg{z) are very similar 
to each other. This is shown in Figure 3 for simulations with the numerical model of 
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Figure 3. Tg/Tgo (solid lines) and p/po (dashed lines), where Po = /o(0), plotted 
versus height for Earth conditions with d = 250prn and two different shear velocities, 
M* = 0.3m/s (blue) and u, ~ O.Sm/s (red). 

Kok and Renno [27]. Since the profile p{z)/M describes the hop height distribution of 
saltating grains, Zs is related to the saltation height h, which is the average hop height 
of the grains. This is discussed in detail in the following section. 

2.7. Physical meaning of Zs 

It was assumed in previous studies that Zs and the saltation height h are proportional 
to each other [17, 19, 20]. This very natural assumption is however not valid for 
the saltation simulated with the numerical model of Kok and Renno [27], since the 
normalized profiles p{z) / po for = 0.3m/ s and u^, = 0.8m/ s in Figure 3 almost coincide 
with each other at small heights. This means, although Zs is larger for = 0.8m/ s, 
the hop height of grains transported close to the surface is almost the same and hence h 
increases weaker with m^, than Zg. It is therefore reasonable to interpret Zg as the height 
of high-energy saltons and the height Zr up to which the profiles p{z) / po coincide as the 
height of low-energy saltons [19], which remains unchanged with u^. In the following 
we explain the reason for this behavior of p{z). 

In our model p{z) decays approximately exponentially, if and only if {vD does not 
vary much with height z. This can be seen from (23), (26), and (32), which allow us to 
write the differential equation, using pg = —p{v1) (analogous to (44)), 

« = = A ~ -PS, (47) 

whose solution is an exponential decrease, if and only if {vD is constant with z. However, 
there can be a huge difference between the value {v1){0) at the soil, which is fixed by the 
splash-entrainment process, and the value of {vD at larger heights, which is proportional 
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to Zgg. Figure 3 shows indeed a strong deviation of p{z) from the exponential shape 
at very small heights within the low energy layer z <^ Zr- In the light of our analysis, 
this deviation corresponds to a strong increase of {v1){z) from {v1){0) towards about 
ZgC] at larger heights. Note that wind tunnel studies [10, 11] did not notice such a 
deviation from the exponential shape in their measurements of p{z). A possible cause is 
that their lowest measurement points, z ~ 20d and z ~ 40(i, respectively, were already 
too high, and they measured only in a region, where {v1){z) was not varying much 
anymore. Note further that a very similar reason, namely that —{vxVz){Q) is fixed by 
the splash-entrainment process, is the probable cause for the very slight deviation of 
Tg{z) from the exponential shape at very small heights (see Figure 1). Finally note that 
mass flux profiles p{z){vx){z), which are often measured in experiments [8, 9, 41, 42, 43], 
decay slightly weaker than p{z), since the particle velocity {vx){z) increases weakly with 
height. 

After explaining p{z), we now calculate h and z^- First, using (41) and (47) as well 
as partial integration, h can be calculated as 

oo oo 

A = J = -i/ .p(.)d. = ^/ = I = H-"^. (48) 



Since Zr does not depend on and since it must be of the same structure as Zg and 
namely proportional to {vl){z)/g at a typical height z, and since ^ = is the only 
height where {v1){z)/g does not depend on due to the splash-entrainment process, 
Zr must write 

Zr oc ^ = (49) 
9 9 

where we used (14). 

Our analysis confirms the picture of Andreotti [19], who hypothesized that one 
can essentially distinguish two species in aeohan saltation, low-energy saltons (reptons) 
slowly moving in small hops, and high-energy saltons moving fast in huge hops. The 
author hypothesized that Zj. is a measure for the height of the focal-region, a region in 
which steady state wind profiles for different shear velocities intersect. 

In this section we showed that the saltation layer can be characterized by three 
heights. The height of low-energy saltatons z^, the saltation height h, and the height of 
high-energy saltons Zs^ which is also a measure for the thickness of the saltation layer. 
In contrast to the first height, the latter two change with V and thus with u^. The 
prediction of V is therefore subject of the following section. 



2.8. Calculation of z* 

The last step towards the calculation of z^ is to derive an expression for V, the last 
undetermined quantity in (46), our relation for Zs- Zg is the main parameter in our final 
relation for z*, which will be of a similar structure as (4), the relation of Raupach [17]. 
For deriving V we use the following strategy. We first approximately calculate U as the 
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wind velocity at an height Zm, which denotes the height of the mean motion, U — u{zm), 
with the mixing length approximation [18]. Then we compute V by 

V^U-Vr, (50) 

where we use that Vr does not change with (see (33)). 

Following the outhned strategy, we calculate U from the mixing length 
approximation [17, 18, 20, 23, 33] as 

(51) 




2 : 

with 

u{zo) = 0. (52) 

In the absence of sand transport, rg{z) = 0, the mixing length approximation yields 
the undisturbed logarithmic velocity profile ((1) with = Zo). In the presence of 
sand transport, rg{z) ^ 0, the velocity profile deviates from the logarithmic shape. 
In Appendix B u{z) as well calculated based on our model assumption 

of an exponentially decreasing grain shear stress profile, (19), and the calculation 
approximately yields 

ln'-^=(l-!^)ln^-G(^) (53) 
Zo \ u^J 1.78^0 \u*J 

and for z > Q.lZs 

„W = !^l|n4 + 4^E.(£), (54) 
K z* 2ku* \ZsJ 

where Ub is the reduced wind shear velocity at the bed, defined by 

Ub^Ua{0), (55) 

u,{z)^uJl-^, (56) 

V Pwui 

and the exponential integral Ei{x) as well as G{x) are defined by 

oo 

E,{x) = / ^dx' = -0.577 -\nx + Y, L , (57) 

J X a. 

G{x) = 1.154(1 + j; In x)(l -a;)^-^^ (58) 

Note that (54) is only an approximation of u{z) for z > O.lZs- In Appendix B one can 
also find an approximation for z < Zs, which however has a much more complicated 
structure. For the coming calculations, we are only interested in the value U — u{zm), 
for whose calculation both approximations perform similarly well, since Zm is between 
0.12;^ and OAzg as we show later. 

In (53) and (54) the shear velocity at the bed Ub remains undetermined. It was 
shown by Duran and Herrmann [20] that Ub must decrease with u^, starting from 
Ub{ut) = Ut, however with a small slope. The reason is that one observes a focal region 
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in aeolian steady state saltation, called Bagnold-focus, below which the wind velocities 
decrease with increasing u^^ [1, 44]. Andreotti [19] argued that this strong decrease of 
wind velocities, the presence of the Bagnold-focus, and consequently the value of Ub are 
mainly caused by grains transported in the low-energy layer. As we explained in Section 
2.7, at very small heights within the low-energy layer z ^ Zr the mass density profile p{z) 
and also the grain shear stress profile Tg{z) deviate from the exponential decrease. We 
however used this exponential decrease of Tg{z) (see (19)), our main model assumption, 
for the computation of the apparent roughness and the wind profile in (53) and (54). 
This does not mean that (53) and (54) are wrong, because the deviation of rg{z) from 
the exponential shape is very small (almost invisible in Figure 1). But it means that 
we cannot use the real value of Ub, which is influenced by the low-energy layer. Instead 
we must use a value of Ub, which corresponds to the extrapolation of the exponential 
shape of Tg{z) above z = z,. to the height 2; = 0. Such a value of Ub would be larger than 
the real value, because p{z) and thus Tg{z) deviate from the exponential shape towards 
higher values within the low-energy layer. We therefore propose that this value is close 
to Ut, 

Ub ~ ut. (59) 

This is supported by simulations using the numerical model of Kok and Renno [27]. 
Figure 4 shows that wind profiles calculated by (54) with the hypothesis (59) are much 
closer to the simulated profiles than those in which the simulated values of Ub were used, 
and this although this hypothesis eliminates the Bagnold-focus. Especially in the region 
which we are interested in, z > O.lZg, (53) and (54) provide an excellent approximation 
of the simulated wind profile, if using (59). Note that (59) is also known as Owen's 
second hypothesis, and it has been used in many previous models e.g. [6, 17, 21, 22, 23]. 
From (53), (54), and (59) we now obtain 

In^^fl-^Vn^-Gf^l (60) 



and 



[/ = -ln^ + ^V^Ei(7), (61) 



where 7 is the third model parameter and defined by 

7=^. (62) 

Zs 

The first terms on the right hand side of (53) and (60) are identical to the relations of 
previous studies [17, 20] (see also (4)). The second term appears, because we did not 
approximate the right hand side of (51) before the integration as done in these studies. 
Our solution is therefore more precise. Note that the height of the mean motion z^n 
should be in leading order proportional to the decay height Zg of the grain shear stress 
profile, because the mean motion is dominated by the motion of high-energy saltons [19], 
and Zg is a measure for the height of high-energy saltons (see Section 2.7). Consequently 
7 is in leading order constant. 




u / (m/s) u / (m/s) 



Figure 4. u plotted versus z/zg for Mars conditions with d ~ 250 fim and two different 
shear velocities, m* = 0.39m/s (a) and = 0.59m/s (b) . The solid lines shows the 
simulated wind profiles [27], the dashed lines show the wind profiles computed by (53) 
and (54) with Ub = Ut (blue) and the simulated values of Ub (red), respectively. 



Having obtained a relation for U, we can calculate V with (50), insert V in (46) 
to obtain Zg, and insert Zs in (60) to obtain z*. The calculation of can therefore be 
summarized as 

ln^=fl-^Vn^?^-Gf-V (63a) 
Zo V l-78'jZo \u^J 



9 



(636) 



I7=-ln^ + ^^^E,(7), (63c) 

C,{Vr)vl = (63d) 

where we used (60) and (62) for (63a), (46), (50), and (62) for (636), and (63c?) and (63c) 
are the same as (33) and (61), respectively. (636) and (63c) can be solved iteratively 
for U and Zm, and (63 o?) can be solved iteratively given a certain drag law Cd{V). This 
is our novel prediction of z* in the most general version. If the impact threshold Ut is 
known, z* can be calculated using {Q3a-Q3d) as function of the model parameters a, (3, 
and 7. It is furthermore possible to compute the mass flux Q as function of the same 
parameters. For this purpose we first compute the mass of transported sand per unit 
soil area M from (34) and (59) as 

M = ^{ul-u^,). (64) 
9 
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Then the mass flux Q becomes 

oo 

g = / p{z)v^{z)dz = MV, (65) 





Q^^{ul-u',){U-Vr), (66) 

where we inserted (64) for M and (50) for V. U and Vr can be computed by {63b-63d). 
At the moment all model equations are functions of the model parameters and Ut- In 
order to close the model, Ut must be calculated as function of the model parameters as 
well. This is done in the following with the help of a closing assumption. 

2. 9. Closing the model - a relation for Ut 

(63a-63(f) are already a novel expression for z^. Like previous expressions in the 
hterature (see (2)- (4) and (6)) it needs the impact threshold Ut as an input parameter. 

We therefore close the model by deriving a relation for Ut in this section. The strategy 
is as follows. We first motivate a simple expression for the particle velocity at the 
threshold, Vt = V{ut). Wc then compute Vr = Ut — Vt, where Ut = U{ut), and 
combine the result with our previous expression for Vr, {Q3d). The resulting equation 
can be rearranged to compute Uf. 

As outlined before, we motivate the following closing expression, 

Vt = 7jUt + Vo = v-ln^ + Vo, (67) 

K Zq 

where z^t = Zm{ut), K = {Pot'^xot + Po\:Vxo.i) / Po is the average particle slip velocity 
(i.e., the particle speed at the surface), where Vxo-\(\.) = '^2:t(4)(0) ^.nd Pot(4.) — Pt(4.)('-')' 
and 7] is the fourth model parameter. (67) means that the difference between average 
particle and slip velocity under threshold conditions is proportional to the average 
wind velocity. This is justified in the following manner. The particle slip velocity 
K is a quantity that like the model parameters a and /3 only depends on the impact- 
entrainment process. In particular Vg is independent of the average wind velocity Ut- 
Prom theoretical and experimental studies it is known that the profile of the average 
horizontal particle velocity Vx{z) starts with Vg at z = and increases with height z 
[10, 11, 27]. The average increase with z mainly depends on the average wind velocity 
Ut- Consequently, the average particle velocity Vt is a function of the average wind 
velocity Ut plus an offset Vo- The simplest possible relation with such a behavior is 
given by (67). Thereby rj describes how efficiently the wind accelerates transported 
grains under threshold conditions. Rearranging (67), rj can be written as 

. = M. (68) 

We calculate the particle slip velocity Vo with the model of Kok [40], explained in 
Appendix A. The result, Vg as function of gd, is plotted in Figure 5. 
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Figure 5. Vo computed with ttie model of Kok [40] (see Appendix A) plotted versus gd. 
Here g = 9.81to/s^ (solid line) and g = 3.71to/s^ (dashed line) are fixed, respectively, 
whereas d is varied. 



Now we can use (67) to calculate Vr- According to (33), Vr does not depend on 
M*. We can therefore write 

Vr = Vriut) = Ut-Vt={l-v)-ln^- Vo, (69) 

Rearranging and using (46), (62), and (67) finally yields an expression for Ut, which can 
be summarized as 

<Vr + Vo) ^ 

(1 - r/) In ' ^ ' 



Ca{Vr)vl = (70c) 



3. Model Validation 



The model is validated by simulation results with the numerical state of the art model 
of Kok and Renno [27] and by several experiments. The simulation results validate (46), 
(50), (63 a), (63 d), (64), and (69), and confirm our statements that the model parameters 
a and /3 are approximately independent of the shear velocity and atmospheric conditions, 
and further indicate that 7 and rj are not varying much as well. In detail we show that all 
model parameters always adopt approximately the same values for whatever conditions 
are simulated. The experiments confirm our expressions (63a-63(i), (66), and (70a- 70c) 
by showing that the same set of model parameters can explain the following experiments: 
the mass fiux data of Creyssels et al. [10], the apparent roughness data of Rasmussen 
et al. [12], and the combined impact threshold data of different studies [24, 25, 26]. 
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Table 1. Conditions simulated with the numerical program of [27]. The first five 
conditions describe an Earth atmosphere with five different particle diameters d. The 
next four conditions describe a Mars atmosphere with four different values of d. The 
remaining nine simulated conditions are imaginary conditions, where one or more of 
the atmospheric parameters were varied. 



3.1. Independence of the model constants of atmosphere and grain properties 

For the verification of the independence of the model parameters we use simulation 
results of the numerical model of Kok and Renno [27] for conditions, which are 
summarized in Table (1). Fluid densities and viscosities fi as well as particle densities 
Ps are varied between Earth and Mars values and particle diameters d are varied between 
100pm and 500pm. The surface roughnesses Zo in the absence of saltation is chosen to 
be Zo = d/30, except in two cases {zo = d/10 and Zg = d/90) which allow us to check, 
whether the predictive performance of our model equations is sensitive to the value of 
Zo- The model of Kok and Renno [27] uses the drag law of Cheng [37], namely 

which we use to compute Vr in (63 (i) and (70c). 

We evaluate (46), (50), (63 a), (63 d), and (64) in order to show the approximate 
independence of the model parameters a, f3, and 7, of atmospheric conditions and grain 
properties. Exemplary for Earth conditions with d = 500pm and Mars conditions with 
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Figure 6. AI calculated using (64) with a as given in Table 1 versus the values of M, 
obtained directly from the simulation for Earth conditions with d — 500/im (blue) and 
Mars conditions with d — 250/im (red) . Each circle belongs to a different shear velocity 
M,. The solid line indicates perfect agreement. The parameter ut, which appears in 
(64), was not calculated by our model, but directly obtained from the simulations (see 
Table 1). 

d = 250/im, Figure 6 shows M calculated using (64) with a as given in Table 1 versus 
the values of M, obtained directly from the simulations. Furthermore, for the same 
conditions. Figure 7 shows Zs and \n{z*/zo) calculated using (46), (63 a), and (63 d) with 
(3 as given in Table 1 versus the values of Zg and ln{z*/zo), obtained directly from the 
simulations. And finally, for the same conditions. Figure 8 shows V calculated using 
(50), (63 c), and (63 d) versus the values of V, obtained directly from the simulations. For 
all plots, each circle belongs to a different shear velocity and the solid line indicates 
perfect agreement. Conditions different from those plotted in Figs. (6-8) show in most 
cases the same good agreement. The values of a, /3, and 7 for all conditions are given 
in Table (1). It shows that a is between 0.9 and 1 for most of the tested conditions. 
A variance of a and thus M of about 10% is however small compared to the degree of 
uncertainty one usually faces in saltation mass(flux) measurements. Furthermore, /3 is 
between 0.12 and 0.135 for all tested conditions, the variance of f3 is therefore even less 
than the variance of a. Furthermore, Table (1) shows that 7 is between 0.23 and 0.34 for 
all of the tested conditions. Therefore we can confirm that a, (3, and to a lesser extend 7 
can indeed be used as approximately constant parameters for saltation simulated by the 
model of Kok and Renno [27] at least within the range of conditions displayed in Table 
(1). We want to emphasize that in particular (46), which is the main contribution of 
our paper, well describes the behavior of the simulated decay heights Zs (see Figure 7). 
Note that the slight disagreement of the Mars simulations from the perfect agreement 
in Figure 8 is probably due to turbulent fluctuations of the wind velocities, which are 
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Figure 7. Zg (a) and ln(z*/zo) (b) calculated using (46), (63a), and {63d) with [3 
as given in Table 1 versus the values of Zg and ln(z*/zo), obtained directly from the 
simulations for Earth conditions with d = 500/iTO (blue) and Mars conditions with 
d = 250 fim (red). Each circle belongs to a different shear velocity u*. The solid line 
indicates perfect agreement. The parameters ut and a are taken from Table (1) and the 
simulated values of V are used in (46), instead of calculating them with the parameter 
7- 
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Figure 8. V calculated using (50), (63c), and (63d) with 7 as given in Table 1 
versus the values of V, obtained directly from the simulations for Earth conditions 
with d ~ 500/im (blue) and Mars conditions with d = 250/im (red). Each circle 
belongs to a different shear velocity m*. The solid line indicates perfect agreement. 
The parameters Ut and a are taken from Table (1). The parameter /3 is not used, 
instead the simulated values of Zs are taken for the computation of V. 
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Figure 9. Plot of Ut — n^^Uth\{zmt/ Zo) versus Vr + Vo- According to (69) Vr + Vq 
and Ut are proportional to each other with a proportionality constant 1 — rj. Each 
circle corresponds to one of the conditions in Table (1). 

considered by the numerical model of Kok and Renno [27], but not by our analytical 
model. These fluctuations are much more important for small gd like Mars conditions 
with d = 250/im than for large gd like Earth conditions with d = 500/im. 

Finally we check (69) by plotting Vr+Vo over Ut = n~^Ut hi{zmt/zo) for all simulated 
conditions. Thereby Vr is calculated by (63 d) with the values of a in Table (1) and 
the values Vo are also given in Table (1). Note that in the simulations Vo does not 
change significantly with as we also stated before from a theoretical point of view. 
Figure 9 shows that the plotted circles, each of them corresponding to one of the 
conditions in Table (1), approximately lie on a straight line through the origin. This 
indicates an approximately universal behavior of rj, because according to (69) Vr + Vo 
and K.~^Ut lia{zmt/zo) are proportional to each other with a proportionality constant l—rj. 
The values of rj are also given in Tab. (1). 

In conclusion, the simulation results of the numerical state of the art model of [27] 
can be very well described by our analytical model indicating only slight variances of 
all four model parameters with changing conditions. Note that the simulated values of 
a, (3, and Vo do not follow the predictions of the model of Kok [40], which were plotted 
in Figs. (2) and (5). 

We also want to emphasize that we entirely failed to fit the numerical data for 
the apparent roughness z*, when using (6), the scaling for Zg proposed by Andreotti 
[19]. Earth and Mars conditions could not be fitted simultaneously by (4) and (6), 
or alternatively (63a) and (6) with a single proportionality constant in (6). Fitting to 
good agreement with the numerical Earth data, led to a disagreement with the Mars 
data by almost two orders of magnitude, even after trying several modifications like for 
instance replacing y/sgd in (6) by Ut- This strongly indicates that (6) does not describe 
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the physics sufficiently well. Andreotti [19] mentioned himself that (6) is nothing but 
a guess: "We have not found any simple explanation of this scaling with i> = fi/pw 
The only indication that we have identified the good parameter is the prefactor of order 
unity." 



3.1.1. Explicit relation for mass flux Q of monodisperse sand, d = 250/im, on Earth 
Many mass flux relations in the literature are given for free field Earth conditions with 
d — 250//m e.g. [23]. In this section, we provide a further explicit prediction of Q for 
these conditions, based on the parameter values in the third row in Table (1). In contrast 
to the parameter values which we obtained from wind tunnel experiments discussed in 
the following section, the values given in Table (1) correspond to free field conditions, 
because the numerical model of Kok and Renno [27] was adjusted to such conditions. 
From a = 0.94 we obtain Vr = 1.55m/s using {63d) and (71). From Uf = 0.196m/s, 
(3 = 0.125, and 7 = 0.33 we further obtain Zmt = 0.0153m using (636) evaluated at 
= Uf. Inserting all values in (66) then yields 

0.94^1-^ " 



-^ln^ + 2.5Gr^ -1.33 1-^ 



+ 1.05 1 - ^ 



,2 



ul J 



(72) 



where we used (63a-63c). This is however not an explicit relation for Q, since Zm 
increases with as well. In order to obtain an explicit relation, we approximate 
Zm ^ Zjnf This is reasonable, since the increase of Q with z^ is only logarithmic. 
Further inserting Zo = d/30 (see Table 1) and using that the first non- vanishing term 
0.01251n2(l— Mj/m^)^ of G{ut/u^) (sec Appendix B) is very small compared to the other 
terms close to Ut, G ~ 0, we finally obtain 



Qg / u 



= 11-4 



17.67^-1.25 fl- U 0.98 fl-^^ ^ 



(73) 



3.2. Experimental validation 



For the validation with experiments we use the drag law (25) of Cheng [37], which we also 
used before. Furthermore, wc need to account for the fact that the surface roughness 
Zo of a quiescent sand bed is a function of the roughness Reynolds number. This is 
in particular important when comparing with experimental impact thresholds, because 
some experiments were made with very small particle diameters d in the aerodynamically 
smooth regime. The whole context is explained in Appendix C including equations, 
which are used to compute Zo- 

In this section we validate our apparent roughness prediction, {63a-63d), with the 
experiments of Rasmussen et al. [12] for five different particle diameters d. The chosen 
data set has the advantage that the scatter in the data is small in comparison to other 
data sets [13, 14]. Furthermore wc use a combination of several data sets [24, 25, 26], 
in order to evaluate our impact threshold prediction, (70a-70c) and the data set of 
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Creyssels et al. [10] for our mass flux prediction (66). The latter choice is motivated 
by the fact that Creyssels et al. [10] used particle tracking methods, which arc more 
accurate than measurements of the mass flux Q with sand traps, which underestimate 
Q by up to 50% [41, 45]. Furthermore the experiments were performed in the same wind 
tunnel as the experiments of Rasmussen et al. [12] and mass flux measurements typically 
vary from wind tunnel to wind tunnel. For instance the two recent measurements of the 
mass flux with particle tracking methods [10, 11] show the same qualitative behavior, 
an approximate scaling of Q with ul—u^, but quite different magnitudes of Q, although 
they used the same sand in their experiments. We want to strongly emphasize that we 
use only one single set of parameters, a, 7, and rj to flt all data sets at the same time 
and that the values of Uf, obtained from our prediction (70a- 70c), are used in (63a-63(i) 
and (66). 

By fltting the model parameters to a = 1.02, /3 = 0.095, 7 = 0.17, and rj = 0.1 
we obtain good to excellent agreement with all data sets. This is shown in Figs. (10- 
12), which present the comparison of {63a-63d) with the data of Rasmussen et al. [12], 
the comparison of (70a- 70c) with the impact threshold data sets [24, 25, 26], and the 
comparison of (66) with the mass flux data of Creyssels et al. [10], respectively. It 
is important to note that the fltted values of the model parameters signiflcantly differ 
from those in Table (1), which were adjusted to match the numerical data obtained from 
simulation with the model of Kok and Renno [27]. The very likely cause is that the 
numerical model of Kok and Renno [27] was adjusted to field data, whereas our model 
is adjusted to wind tunnel data. According to Sherman and Farrcll [14], field and wind 
tunnel data of the apparent roughness z* can differ by up to one order of magnitude. 
The reason for this is not fully understood. It was speculated by Raupach [17] that 
differences between wind tunnel and field data are attributed to not fully equilibrated 
saltation in wind tunnels. It is however likely that this is not the only reason, since 
even in large wind tunnels like the one used by Creyssels et al. [10] and Rasmussen 
et al. [12], for which we adjusted our model, the differences still exist. We propose 
that differences in the grain size distribution could be another cause. Sand in field 
experiments has typically a much broader grain size distribution than sand with the 
same mean diameter d typically used in wind tunnels. For instance the five different 
sands used by Rasmussen et al. [12] are the same as the ones used by Iversen and 
Rasmussen [46], which the latter described as "closely sized sand samples" and referred 
to as "uniform samples". Bagnold [21] argued that the average particle velocity V in 
aeolian saltation is found to be about two times higher for broadly in comparison to 
narrowly distributed grain sizes with the same d, because "elastic rebound of smaller 
grains off more massive ones is superimposed on the splash process", leading to higher 
saltation heights. According to (46), two times higher V in field experiments than in 
wind tunnel experiments would lead to about three times higher Zg and therefore to a 
much higher apparent roughness z^. 
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Figure 10. Comparison of our model (63a-63(i) (solid lines) with experimental data 
of Rasmussen et al. [12] (circles) for five different particle sizes, (a) d = 125/^to, (b) 
d = 170/im, (c) d = 2A2fj,m, (d) d = 320/zm, and (e) d = 544/ito. (c) includes six z* 
data points (red circles) measured by Creyssels et al. [10]. 
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Figure 11. Comparison of our model (70a-70c) (solid line) with experimental data 
[24, 25, 26] (circles). 

0.1 1 ' ' ' ' 1 
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Figure 12. Comparison of our model (66) (solid line) with experimental data of 
Creyssels et al. [10] (circles). 

4. Discussion 

We have presented a comprehensive analytical saltation model, which provides a new 
level of understanding of the change of the average wind profile during aeolian saltation 
and of aeolian saltation in general. The model incorporates (63a-63rf), a novel expression 
for the apparent roughness z*, and (46), a new expression for the thickness of the 
saltation layer Zg. (46) can be seen as the main contribution of our study and it is to 
our knowledge the first relation for Zg entirely derived from physical principles. It is 
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an improvement of the relation (6), proposed by Andreotti [19], which had no physical 
foundations. Besides, our saltation model also provides (66), a novel expression for the 
mass flux Q, and (70a-70c), an expression for the impact threshold Uf. All relations 
are extensively evaluated. They are in very good agreement with the recent numerical 
model of Kok and Renno [27] and with many experiments [10, 12, 24, 25, 26]. 

In particular our expression for the mass flux Q, (66), is in excellent agreement 
with the data set of Creyssels et al. [10] (see Figure 12). The authors found a scaling of 
Q with ul — u^, and it can be seen that this is also the dominant term in (66). The same 
scaling was also found by theoretical studies [19, 34, 47, 48], and recently measured by 
Ho et al. [11], who showed that this scaling relation is a consequence of the bed being 
erodible instead of fixed. The splash-entrainment mechanism on erodible beds causes 
the particle slip velocity Vo to be constant with increasing m*, what in turn hinders the 
average particle velocity V to increase fast with u^,. The main increase of the mass flux 
Q = MV with u^, comes instead from the average transported mass per unit soil area 
M, which scales with ul — u^, as we and other studies e.g. [23] showed. 

In conclusion, the facts that our model stays on sound physical foundations and 
that it is in very good agreement with experiments and state of the art simulations 
make us confldent that it for the flrst time rehably quantifles the feedback of the sand 
transport on the wind momentum, even on planets different from Earth. 
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Appendix A. Computing impact and lift-off velocities with the model of 
Kok [40] 

Kok [40] has recently derived analytic expressions for the average impact and lift-off 
velocities, from which wc obtain the following expressions for the constants a' and 
and for the particle slip velocity Vg. They write 

, -Avzo sin diVi + sin 9iVi , . ^ . 

^ = ~ii = 2 — > l^-lJ 

/\Vxo cos UiVi — cos OiVi 

^ Avl {cos OiViY - (cos eiviY' ^ • ' 

_ Po^Vxot + PoiVxoi _ sin 9iVi cos 9iVi + sin 9iVi cos OiVj 
Pot + Pol sm OiVi + sm OiVi 

where 9i ^ 11°, 9i ^ 40°, and Pot'^^zot + Poi^zoi — were used (see (14)). Vi and vi are 
furthermore given by 



1-F , — - 1 11 fi-rV , i + F 



^^-^^9fd-^^ + dj^, + [^ 9fd+^V9fd, (A.4) 
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and 

VI = Fanv.^, (l - (j:^) + "e.-. (^1 - exp (-^^J 1 , (A.5) 

where the parameters are given hy r — 0.02, F — 0.96, e ^ Is/m, — 0.55, 
a^j = Q.lb^g/gf and 5^/ is an effective gravity, which incorporates the effect of cohesion 
at small particle diameters 0?, 

»' = »+^- 

Here C, is the dimensional cohesion parameter. (A. 6) expresses that cohesive forces scale 
with in contrast to the gravity force, which scales with d^, which means that at small 
d cohesive forces become dominant. Shao and Lu [49] estimated the magnitude of C, to 
be between about lxlO~^A^/m and 5xlO~^A^/m, we use C, — 5xlO~^A^/m. 

Appendix B. Calculation of the average wind velocity profile and the 
appeirent roughness 

In this appendix we show how one can compute the average wind velocity profile and 
the apparent roughness from an exponentially decaying grain shear stress profile, 

r,(^) = ve-^/^^ (B.l) 

The average wind velocity profile u{z) can be obtained from Prandtl's mixing length 
approximation [17, 18, 20, 23, 33] as 

d.u{z) Ua{z) 



with 

u{zo) = 0, (B.3) 

where r = pwul- The apparent roughness 2;*, is the roughness of the asymptotic fluid 
velocity profile v{z) = v{z)\z^oo, giving 

M = — In-. (B.4) 

In order to integrate (B.2) analj^ically, we Taylor-expand the square root in (B.2) in 

the argument aexp{—z/zs), where a — Tgo/r. It becomes 

00 

VT^^l-Y,fP'^ (B.5) 

where fj is given by 

(2j - 3)!! 

S> = (B.6) 
Then (B.2) becomes 
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Within the integral we transform the z-coordinate using 

x = ^—, (B.8) 
such that the integral transforms to 

Z jz/Zs 

J^^-jz'/z.^ J ^e-.^Ei(jV^,)-Ei(jV^,), (B.9) 

Zo jZo/Zs 

where Ei{x) is called the exponential integral. It can be expressed as 

Ei(a;) = -0.577 -In a; + Ein(x), (B.IO) 
where 0.577 is the Euler-Maschcroni constant, 

0.57721... = lim ^ - ln(n) j , (B.ll) 

and Ein(a;) is an analytic function with the series-expansion 

1=1 

Ei{x) vanishes for a; ^ oo. Using (B.9) the velocity profile finally writes 



oo 



uiz) = ^[ln--J2 fja' (Mj^o/zs) - E,ijz/zs)) . (B.13) 
The comparison with the asymptotic profile (B.4) then yields 



oo 



ln^ = J]/yEi(jV^,). (B.14) 

Since Zg <^ Zs, (B.14) can further be written as 

In ^ = -K{a) + (1 - ^r^) (in - - 0.577 ) , (B.15) 

\ Zq J 

where K[a) is defined by 

oo 

K{a) = Y,fM3y- (B.16) 

To our knowledge, the sum on the right hand side is not analytically treatable with 
common methods. However, we found that it can be very well approximated by 

K{a) ^ 1.154(1 - v^r^)2-56(i + Vr^ln^r^), (B.17) 

where both constants, 2.56 and 1.154, are fit constants, ensuring very good agreement 
between the infinite sum and the approximation. This approximation performs very 
well over the whole range of a, the relative errors being typically below 1%, as is shown 
in Figure Bl. Now we can finally write 

ta£^=fi_!^!;)t^ + Gf^V (B.18) 

where Ub — ~ « was used (see (55)), 1.78 = exp(0.577), and G{x) is defined by 

G{x) = 1.154(1 + X In a;)(l -a;)^-^^ (B.19) 




Figure Bl. Comparison between K{a) computed using (B.16) with an accuracy 
of at least 10^^ (squares) and K{a) computed witli tlie approximation (B.17) for 
< a < 0.999. 



Appendix B.l. Approximation for large z 

Now we further calculate an approximative expression for the wind profile at large 
heights z. The velocity profile in (B.l 3) can be rewritten as 

u{z) = ^ l^ln^ + fj/,a^Ei(j^/z,)j . (B.20) 

where we inserted (B.14). If z/zg is large enough, it is sufficient to only consider the 
first term of the Taylor-expansion, because Ei(a;) decreases strongly with increasing 
argument x. This eventually yields 

«W = ^^llu4 + 4^E,(±). (B.21) 

where we also inserted (55) and (B.18). 

Appendix B.2. Approximation for z < Zg 

Here we motivate an approximation of the wind profile u{z) for small heights z < z^. 
For this purpose, we first make an approximation for the infinite sum 

(\ °° 
^'^j :=5Z/yEin(jV^,), (B.22) 
i=i 

and then insert it into the velocity profile, given by (B.20). Inserting the series-expansion 
(B.12) of Ein in (B.22) and exchanging the order of the sums gives 

k=l j = l 
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The inner sum can be written as 

oo 

It can be shown that 

k 

1=1 

where the Stirhng numbers of the second kind Sj^^ are defined by 

^f = ^E(-irQ('-™)'- 

m=0 ^ ^ 

Using 

Sill - ^] = 

we can write 

^ - ^ W. ^ 2'(l-ay-V2- (^-24) 

k=l 1=1 ^ ' 

The above form of / and the fact that / can be approximated by the first-order Taylor- 
expansion 1 — |aEin(^/^s) for small a allows the following approximation, 

m 2^(l-a)'=-V2' ^""-^^^ 

which is the Taylor-expansion of the function 

aEin(^A,) 3F2 (i,l,l;2,2,-^5ig^) 

2v^r^ ^ ' 

The function 3F2 is a generahzed hypergeometric series [50]. This approximation 
performs very well for z < Zg and the whole range of a, as can be seen in Figure 
B2. The relative errors of this approximation are typically below 1%. Even for larger 
values of z up to 102;^ the approximation is still good, with relative errors below 10%. 
We can rewrite (B.26) in more compact form as 

I = ( ~ ul)^ in{z/zs 

2m* V ^6 
where we used (55) and the function H[x) is defined by 

//(x) = X3F2(^,2,2;2,2,-x). (B.28) 

If accuracy is not crucial, the complicated function H[x) can be replaced by x^''^^ for 
small arguments x < 20. This is shown in Figure B3. As outlined before, we now insert 



-H ( ^"^"-^ ) , (B.27) 
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Figure B2. Comparison between / as function of z/zg computed using (B.22) with an 
accuracy of at least 10^^ (squares) and / computed using the approximation (B.26). 
The different lines correspond to a = 0.1 (blue), a = 0.5 (red), a — 0.8 (green), and 
a = 0.999 (brown). 
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X 

Figure B3. The function H{x) (solid line) compared with the function (dashed 
line) . 

the approximation (B.26) in (B.20). By further using Zo <C Zg, the approximated wind 
profile u{z) can be written as 

The approximated wind profile calculated using (B.29) (solid line) and (B.21) (dashed 
line) is exemplary plotted in Figure B4 versus the exact solution of the initial boundary 




u / (m/s) 



Figure B4. The exact solution u{z) of tlie boundary value problem (squares), 
calculated using (B.13) with an accuracy of at least 10~^, plotted versus the 
approximation (B.29) (solid line) and approximation (B.21) (dashed line). Here 
d = 250/Ltm, — 0.7m/ s, Ub — 0.1m/ s, and — lOOd are used. 

value problem (squares), calculated by (B.13) with an accuracy of at least 10^^. It can 
be seen that (B.29) is an excellent approximation of the exact solution of the boundary 
value problem for small z and (B.21) an excellent approximation for large z. Since 
both approximations underestimate the analytic solution, the maximum value of u{z) 
calculated using (B.21) and (B.29), represents an excellent approximation for the whole 
range of z. 

Appendix C. Surface roughness of a quiescent sand bed 

The surface roughness Zo in the absence of saltation depends on the roughness Reynolds 
number 

where ks is the equivalent Nikuradse roughness [52, 53]. kg equals the grain diameter d, if 
the grains of the sand bed are monodisperse, spherical, and very well arranged, meaning 
that the center point of each particle of the topmost layer is at the same height. However, 
under more natural conditions kg can be larger, depending on the grain size distribution 
and the arrangement of the sand bed, and the shape of the grains. A typical value, 
which is used by engineers, is kg = Sdgi for water flows in pipes and flumes [53], where 
dsA denotes the diameter value which is larger than 84% of the grains of the grain size 
distribution, or kg = 2d for wind flows [15]. However, since we validate our model with 
experiments, which used narrowly distributed sand [10, 12], we use kg = d. Note that 
the value of ks does not much influence the final results in most cases. The well known 
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Figure CI. Plot of Zo/kg over R^w according to ttie relation of Cheng and Chiew 
[51], (C.4) and (C.5). 



and widely used roughness law 

Zo = K/30 (C.2) 

is obtained for large roughness Reynolds numbers Rew > 100, which is called the 
aerodynamic rough regime. On the other hand, in the limit of low roughness Reynolds 
numbers Rew < 3, the roughness is proportional to the size of the viscous sublayer 

Zo = fi/{9pwU^), (C.3) 

which is called the aerodynamic smooth regime. Most natural conditions for aeolian 
saltation on Earth fall between those regimes, what is called the aerodynamic 
transitional regime. For instance we obtain Rew ~ 6 for wind with a shear velocity 
of = OAm/s over a typical sand surface with a mean diameter of ci = 250 fim and 
kg = d. The behavior of Zq as function of Rew was measured by Nikuradse [52] for pipe 
flows and described by Cheng and Chiew [51] by the following empirical relation 

= /Cs exp(— kI?), (C.4) 

where 

B = 8.5 + (2.51ni?e^ -3)exp {-0.11(ln Rewf '^) ■ (C.5) 

Zo/kg is plotted in Figure CI as function of Rew- It describes a function, which has a 
minimum in the transitional regime at Rew = 9.6 and converges against (C.2) for large 
Rew and against (C.3) for low Rew- The same behavior was measured by Dong et al. 
[54] in wind tunnels. Since the variance of Zo/kg between the transitional and rough 
regime is not very large and the measurement errors of Zo are large, it is appropriate 
for saltation models to use a constant value for Zo/kg like Zo = kg/ 30 in these regimes. 
However one cannot neglect the very strong increase of Zo for Reynolds numbers below 
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Rew = 3 in the smooth regime. Since we consider very small particle diameters in our 
model, for which R^yj < 3, we use (C.4) and (C.5) to compute Zg- 
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